/*  sdflogxe.cpp  by K.Tsuru  */
// function ID 3404     DRADIX  Renewal since ver. 2.30
/********************************************************************
SDouble class
It returns the natural logalithm log x for x > 0 using Exp(x).
[Algorithm]
It is well known that the convergence of logarithmic series is very slow. It costs
much time in a normal method. Then we shall introduce the following idea.
Pay attention to an identity
     x = exp(r){1-(1-x*exp(-r))}
we obtain
     log(x) = r+log(1-x1)
where x1= 1-x*exp(-r).
The value of 'r' can be arbitrarily chosen because of identity.
The solution of an equation x1=0 is r=log(x), then evaluating it by double
the 'x1' has a very small value, order of 10^(-15).
Further solving a equation (1+u)/(1-u)=1-x1 with respect to 'u' we obtain
    u = -x1/(2.0-x1);
and
    log(1-x1) = log{(1+u)/(1-u)}= 2*(u + u^3/3 + u^5/5 +...).
The square of 'u' has an order of 10^(-30) for any x(or r), then in order to
obtain the value in 400 figures it is necessary to sum up only about 15 terms.
In the practical program an intermediate step is inserted in which the approximate
value of 'r' is evaluated in rather large figures.The processing speed depends
on that of exponential function Exp(x).When the value of x is too large or too
small to treat by double, the processing time of Log(10) is added if it has not
been evaluated yet.In othe case it does not depend on the value of 'x'.
*********************************************************************/
#ifndef SN_H
#include "sn.h"
#endif
static void LogApproX(const SDouble& x, SDouble& approX);   //sub-function
static const char* const func = "LogE";
static const int appFig = 32;
// static const SDouble ONE(1.0);

SDouble LogE(const SDouble& x){
    SDouble X, add;

// x = X*10^exp. log(x) = log(X)+exp*log(10), add = exp*log(10), 0< X < 1
    if( GetLogxCalcMethod(x, X, add) ) return X; // x= 1 or 10, X=0 or log(10)

    RealSize C;

// step 1 :It gets an approximate value by use of double log(double).
    uint ef = x.EffFig();
    uint fig = min( (uint)appFig, ef );

    //If the value of 'fig' is increased with figures, the speed is almost the same.
    C.SetEffFig(fig);
    SDouble r;

    r.SetZero();
    LogApproX(X, r);
    C.SetEffFig(0);

// step 2 :Using the approximate value of step 1 it gets in the full precision.
    if(ef > (uint)appFig) LogApproX(X, r);
    if(add.Sign(3401)) r += add;

    return r;
}
/*********************************************************
Zero or approximate value is given to 'approX'.
When approX=0 it takes as r=log(doubleD(x)). In other case
r=approX.
The value
log(x) = r+log(1-x1)
evaluated by series is overwritten on 'approX'.
*********************************************************/
static void LogApproX(const SDouble& x, SDouble& approX){
    SDouble x1, r, delta;
    delta = ONE - x;

// x == 1.0 ?
    if( delta.Sign(3401) == 0 ){    // X == 1.0
        approX.SetZero();   return;
    }

    int fexp;   //fixed exponent when r != 0, r = log(x)
    if(delta.RdxExp() < -appFig){
//When x is very close to 1 (x==1 in double precision), xd = 1.0, r = 0.0 and x1 = 1.0 - x.
//It includes the case in which 'doubleD' cannot be used.
        r.SetZero();    // r = 0
        x1 = delta;     // delta = one - x;
        fexp = approX.Sign() ? approX.RdxExp() : x1.RdxExp();
    } else {
        RealSize C;
        uint upPrec = 0;
        if(approX.Sign(3401)) r = approX;
        else{
            r = log(doubleD(x)); //evaluated by double
/*
Because |x1| is very small it increases the effective figures to avoid
the cancellation. This is necessary to exactly obtain the self-evident
relation "log(e^1.25) = 1.25" by the following judgement "if(x1.IsMLT(r))".
*/
            upPrec = x.Hidden();
            upPrec = x.ProperUpPrec(upPrec);
            if(upPrec) C.SetEffFig(x.EffFig() + upPrec, C.TEMP_EXTEND);
        }
        //It evaluates by SDouble, the sign of 'x1' is unknown.
        x1 = ONE - x*Exp(-r);
    // |x1| is O(10^(-15))
        if(upPrec) C.SetEffFig(0);
        if(x1.IsMLT(r)){    // 'x1' is a value including error only.|x1|<<|r|
            approX = r;     return; //A precision value has been obtained because
                                    //of the nice value of 'approX'.
        }
        fexp = r.RdxExp();
        if(x1.RdxExp() > -2) x.SetError(x.FATAL, func, 3401);
    }

/*
Evaluation log(1-x1) by series.
It enters into the fixed point mode with a fixed exponent 'fexp'.
Here the vlaue of 'x' is short or very small decimal fraction, then
it does not provide the SDecimal function.
See also "SinSeries" and "Asin".
*/
 // (1+u)/(1-u) = 1 - x1, x1 and u is very small.
    const SDouble u = x1/(x1 - 2.0), usq = u*u;
    r.FixedPoint(fexp);

#if 1 // 5.85 (sec)
    SDouble v = u, sum = u; // sum = (1/2)log{(1+u)/(1-u)}= u + u^3/3 + u^5/5 +u^7/7+...
    v *= usq;

    sum +=  DsDiv(v, 3);
    ulong k = 5, mt = x.SlOpMaxValue();

    do{
        v *= usq;
        delta = DsDiv(v, k);
        sum += delta;
        if( (k += 2) > mt ){
            sum.SetError(sum.NOT_CONVERGE, func, 3401);
            break;
        }
    } while(delta.Sign(3401));
    r.PointFree();
    approX = DsMult(sum, 2);    // includes Reform()

#else  // 5.87(sec)
    SDouble sum = ONE;      // log{(1+u)/(1-u)}= 2u*(1 + u^2/3 + u^4/5 +u^6/7 +...)
    SDouble v = usq;
    sum +=  DsDiv(v, 3);
    ulong k = 5, mt = x.SlOpMaxValue();

    do{
        v *= usq;
        delta = DsDiv(v, k);
        sum += delta;
        if( (k += 2) > mt ){
            sum.SetError(sum.NOT_CONVERGE, func, 3401);
            break;
        }
    } while(delta.Sign(3401));

    r.PointFree();
    approX = u * DsMult(sum, 2);    // includes Reform()
//end of the calculation of 'log(1-x1)'
#endif

    if(r.Sign(3401)) approX += r;
    x.upToTerm = (k-1)/2;
    return;
}